%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function computes the discrete state markov chain approximation for
% a continuous state uni-variate AR(1) process using Tauchen (1986) method.
% Function: 
% [StateGrid,TrProbMat,InvDist]=TauchenScalar(GridWidth,GridNum,rho,sigma,TrProbEps)
% Input variables: 
%   GridWidth:  Width of the grid (in the unit of ergodic dist. std.)
%   GridNum:    Number of the discrete state space grid
%   rho:        AR(1) coef. of the process
%   sigma:      std. of the AR(1) process innovation term
% Output variables:
%   StateGrid:      discrete state space grid
%   TransProbMat:   transition prob. matrix
% Special Note: TransProbMat_{i,j}=Pr(Z1=Zi,Z0=Zj)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [StateGrid,TrProbMat,InvDist]=AR1Approx_Tauchen(GridWidth,GridNum,rho,sigma,TrProbEps)

%% Test Part
% rho         =   0.9;
% sigma       =   1;
% GridNum     =   9;
% Width       =   6;
if nargin<=4
    TrProbEps   =   1e-6;
end
%% Computation Part

% std. of the Ergodic Distribution
ErgoSigma   =   sqrt( sigma^2/(1-rho^2) );
% Grid for State Space
StateGrid   =   linspace(-1,1,GridNum)'*GridWidth/2*ErgoSigma;
% Half of the Grid Interval Distance
dZ          =   GridWidth*ErgoSigma/(GridNum-1)/2;
% Construct the Transition Prob. Matrix
Zi          =   repmat(StateGrid,[1,GridNum]);
Zj          =   repmat(StateGrid',[GridNum,1]);

TrProbMat   =   normcdf(Zi+dZ-rho*Zj,0,sigma)-normcdf(Zi-dZ-rho*Zj,0,sigma);
TrProbMat(1,:)    =   normcdf(Zi(1,:)+dZ-rho*Zj(1,:),0,sigma);
TrProbMat(end,:)  =   1-normcdf(Zi(end,:)-dZ-rho*Zj(end,:),0,sigma);
TrProbMat(TrProbMat<TrProbEps)    =   0;
TrProbMat   =   TrProbMat./repmat(sum(TrProbMat,1),[GridNum,1]);
TrProbMat   =   sparse(TrProbMat);
% Construct the Invariant Distribution
[InvDist,~] =   eigs(TrProbMat,1,1+1e-6);
if sum(InvDist)<0
    InvDist     =   -InvDist;
end
InvDist(InvDist<0)      =   0;
InvDist         =   InvDist/sum(InvDist);

